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Abstract 

An analytical procedure is presented, called the 
modal element method, that combines numerical grid 
based algorithms with eigenfunction expansions devel- 
oped by separation of variables. A modal element 
method is presented for solving potential flow in a 
channel with two-dimensional cylindrical like obstacles. 
The infinite computational region is divided into three 
subdomains; the bounded finite element domain, which 
is characterized by the cylindrical obstacle and the 
surrounding unbounded uniform channel entrance and 
exit domains. The velocity potential is represented 
approximately in the grid based domain by a finite ele- 
ment solution and is represented analytically by an 
eigenfunction expansion in the uniform semi-infinite 
entrance and exit domains. The calculated flow fields are 
in excellent agreement with exact analytical solutions. 
By eliminating the grid surrounding the obstacle, the 
modal element method reduces the numerical grid size, 
employs a more precise far field boundary condition, as 
well as giving theoretical insight to the interaction of the 
obstacle with the mean flow. Although the analysis 
focuses on a specific geometry, the formulation is general 
and can be applied to a variety of problems as seen by 
a comparison to companion theories in aeroacoustics and 
electromagnetics . 

Introduction 

Computational fluid dynamics (CFD) now plays a 
major role in aeronautical research and design. 
Pironneau (1989) reports that for Dassault industries 
1986 was the year when the numerical budget over took 
the budget for experimentation in wind tunnels. This 
comparison of budgets is appropriate since both wind 
tunnels and CFD analysis provide similar information 
about the flow physics. Numerical grid based solutions 
with color graphics and computer generated movies are 
closely akin to their experimental counterparts with 
smoke streamlines, shadowgraph, and schlieren photo- 
graphy. Current CFD programs in aeronautics model 
with incredible accuracy the flow field of whole aircraft 
as well as propulsion systems and rotor dynamics 
(Henne, 1990). Nevertheless, current CFD analysis does 
not develop the mathematical insight into the role of 
key variables and parameters that readily unfold from 


exact and approximate analytical solutions. Lighthill 
(1989) addressed this inherent deficiency when he stated 
that it is “essential to stress the connections between 
theoretical analysis on one hand and experimental 
and/or CFD studies on the other.” 

To gain more theoretical insight in CFD problems, 
the present paper presents an analytical procedure, 
called the modal element method, that combines numer- 
ical grid based algorithms with eigenfunction expansions 
derived by separation of variables. Herein, a modal ele- 
ment method is developed for solving potential flow in 
an infinitely long channel with a two-dimensional cylin- 
drical like obstacle present. Emphasis is on the problem 
formulation. Wholly numerical finite element solutions 
for the streamlines and potential lines over a cylinder in 
a duct are presented in many introductory texts on 
finite elements, such as Hinton and Owen (1779, p. 222) 
and Segerlind (1976, p. 183), as well as advanced fluid 
dynamic texts such as Chung (1978, p. 177). However, 
the modal element method adds theoretical insight to 
both the numerical formulation and the physical 
problem. In the numerical analysis, the method will aid 
judgment in choosing the grid density as well as the 
accuracy of the exit boundary condition. In the flow 
problem, the method will determine the physical 
parameters which dictate the change of the flow 
streamlines and the potential lines. 

Although the analysis focuses on a specific geometry, 
the formulation is quite general and can be applied to a 
variety of problems as seen by a comparison to its 
companion theories in acoustics and electromagnetics. In 
CFD applications, however, singularity requirements 
introduce some differences from the previous wave prop- 
agation formulations. Historically, a primary reason for 
developing the modal element method was to accurately 
describe the radiation boundary condition at the compu- 
tational boundary of a numerical grid. In electro- 
magnetics, Chang and Mei (1976), Lee and Cendes 
(1987), and Baumeister and Kreider (1993) applied the 
method to scattering from dielectric cylinders while 
Baumeister (1991) applied the method to electro- 
magnetic propagation in ducts with surface irregular- 
ities. In acoustics, Astley and Eversman (1981) 
employed the method in duct propagation problems 
while Baumeister and Kreider (1992) have applied the 
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technique to scattering from soft and rigid bodies* For 
validation, numerical calculations using the modal ele- 
ment method for sound propagation in a variable area 
duct with a cylindrical like obstacle show excellent 
agreement with experimental results (Baumeister, et al., 
1983). 

To illustrate the advantage of coupling analytic and 
grid based numerical solutions in the modal element 
method, consider the problem of finding the pressure 
amplitude resulting from scattering of an acoustic plane 
wave by a rigid cylinder. As shown in Fig. 1(a), a 
conventional finite difference theory (Khan, Brown, and 
Ahuja, 1986) requires a large dense grid to resolve the 
wave like nature of the pressure field and to accurately 
approximate the far field radiation boundary condition. 
In contrast, the limiting case of modal element method 
for rigid bodies requires only a single line of elements as 
shown in Fig. 1(b). Figure 1(c) compares the pressure 
amplitude of these analyses with the exact theoretical 
results shown by the solid line. The modal element 
method shows excellent agreement with the theory while 
the conventional finite difference theory shows some 
error because of the previously alluded approximations. 
The modal element method can also be extended very 
easily to higher frequencies, as shown in Fig. 1(d), since 
no nodes are required in the far field. 

Considering open systems as in Fig. 1, the modal 
element method is a grid based numerical system that 
has many advantages of the classical boundary integral 
methods such as the boundary element method in acous- 
tics, the panel method in aerodynamics, and the method 
of moments in electromagnetics. These boundary inte- 
gral methods are well suited for solving the scattering 
problem discussed in Fig. 1. However, for the semi- 
infinite duct problem considered herein, the boundary 
element method requires a closure approximation in the 
far field similar to the standard finite element method 
(Brebbia, 1978, p. 80), Abo, similar to the modal 
element method considered herein, the finite element 
method and the boundary element method can be com- 
bined as discussed by Brebbia (1978, p. 178) but with- 
out the modal element advantage of obtaining a closed 
form analytical solution. Yet in a broader sense, the 
modal element method could be considered a subset of 
the boundary element method under the title of indirect 
method of analysis (Brebbia, 1978, p. 2). 

The motivation for adapting the modal element 
method for CFD analysis herein is threefold: first, to 
explicitly determine the importance of physical param- 
eters by obtaining a closed form analytical solution in 
part of the solution domain, second, to minimize the size 
of the regions requiring numerical grids, and thirdly, to 
more accurately approximate the exit boundary condi- 
tion of the numerical grid. The later consideration can 


be important in the direct computation of aerodynamic 
sound from unsteady flows where the sound levels can be 
three orders of magnitude smaller than mean flow 
variables. 

In the present paper, first the method of analysis and 
domain decomposition is discussed followed by a devel- 
opment of the analytical solution. Next, subdomain 
interface conditions are presented followed by the finite 
element solution procedure and the requirements of a 
Dirichlet boundary condition. The geometric model and 
an exact solution for the model are presented next fol- 
lowed by results and comparisons that include two 
example calculations. 

Nomenclature 

A total dimensionless area of finite element domain 
A^ area of element e 

A” modal amplitude, Eq. (14) 

B* modal amplitude, Eq. (17) 

b dipole strength, Eq. (37) 

c” modal amplitude, Eq. (9) 

c+ modal amplitude, Eq. (9) 

F column vector, right hand side of global matrix 
equation 

h dimensionless channel height, h^/L^ 

i yTl 

K global matrix, Eq. (33) 

k# wave number, Fig. 1 

characteristic distance 
m mode number 

m* mode number, Eq. (19) 

M number of elements in finite element domain 
M coe f number of modal coefficients used in 
eigenfunction expansion 

M pts number of grid points on interface S used in 
integration, Eq. (20) 

N number of nodes in finite element domain 

Ng number of nodes on S to resolve harmonics 
N^ local linear triangular interpolation function, 

N (e) (x,y); N| e )(xj,yj) = (i = 1,2,3; j = 1,2,3) 
n outward unit normal vector 

n normal vector 
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Rj global residual error at node I 
r^ radius, Fig. 1 

S line interface about finite element domain 

S a entrance plane 

S b exit plane 

S + region exterior to S 

S — region interior to S 

s arc length parameter on S 

U* free stream velocity 

normalizing velocity 

Wj global weight function associated with node i; 

W I (x, ) y J ) = f IJ (I = 1...N; J = X...N) 
wj e ) local weight function associated with node I 
W m . interface weight function, Eq. (19) 

W(z) complex potential, Eq. (37) 

X separation function, Eq. (4) 

x dimensionless axial distance, x^/L^ 

x in starting position of finite element grid 
x Q axial intercept of obstacle 

x out position of finite element grid 

Y separation function, Eq. (4) 

y dimensionless transverse distance, y^/L^ 

y 0 height of obstacle 

z complex variable, Eq. (38) 

6 angle between element outward normal and 

positive x axis 
Laplace operator 

Kronecker delta (<$jj == 1 for I = J; £jj = 0 for 
I * J) 

A eigenvalue, Eq. (8) 

$ column solution vector, Eq. (33) 

<f> potential 


< f> potential 

ip stream function 

Subscripts 

a analytical solution 

b analytical solution 

I global node index in finite element domain 
I s global node index for interface S 

i local element mode index 


J global node index, Eq. (25) 

j local element node index 

Superscripts 

approximate solution 
# dimensional quantity 

( ) average value 

(e) element value 

Method of Analysis 

The present study is concerned with computing the 
potential flow field in a channel with two-dimensional 
cylindrical like obstacles, as shown in Fig. 2. A uniform 
velocity profile U* is assumed to exist far upstream. For 
inviscid and irrotational two-dimensional flow, the 
potential equation governs; 

v* 2 <t>* = o W 

where # denotes a dimensional variable. For this paper 
the following dimensionless variables are introduced: 


X# y 


L* " u 


U # L* 

O 

The superscript -f indicates the direction of the velocity 
in the -h x direction. All symbols are defined in the 
nomenclature. Equation (l) becomes 

( 3 ) 

V 2 <f> = 0 

The exact shape of the obstacle is defined by an 
infinite row of doublets transverse to a uniform flow 
(Kirchhoff, 1985). For obstacles less than half the height 
of the duct, its shape is nearly a circular cylinder. The 
advantage of using this obstacle is that an exact solution 
exists for validating the theoretical results. The detailed 
shape of the obstacle will be full described in a later 
section of this report. 

The common numerical approach to this problem is 
to extend the grid far upstream and downstream of the 
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obstacle such that the assumption of uniform velocity is 
valid, as shown in Fig. 3. Generally, a large entrance 
and exit grid region must be selected. The leading edge 
of the obstacle is defined by x Q while the inlet and exit 
positions of the grid are defined by x^ and x out 
respectively. Symmetry about the cylinder was not used 
allowing the computer program to handle a greater 
variety of problems. 


v = - = 0 at y ~ 0 and y = h (®) 

dy 

yields 

= c m cos(Ay)X(x) ( 7 ) 


In contrast to the conventional approach, a typical 
modal element grid system is shown in Fig. 4. The 
spatial domain is divided into three subdomains, the 
entrance and exit analytical domains and the finite 
element domain. The finite element domain contains a 
nodal grid system that covers the region of complicated 
geometry. Linear triangular elements are used and the 
subdomain interface is approximated by piecewise linear 
segments. In the finite element domain, an approximate 
solution for the velocity potential is calculated by the 
Galerkin method. In the analytical domains, which 
extends to ±«, an eigenfunction expansion for the 
velocity potential is derived by separation of variables. 

The modal element method couples the analytical 
and numerical solutions by imposing continuity on the 
potential and velocity at the interface between the 
subdomains. This coupling results in a single matrix 
equation in which the eigenfunction coefficients and the 
potential at the finite element nodes are calculated 
simultaneously, yielding a global representation of the 
potential field. Next, Eq. (3) will be solved by the 
method of separation of variables in the entrance and 
exit regions and then by the finite element technique in 
the grid system surrounding the obstacle. 

Analytical Solution 


where the eigenvalues are 

A =— m = 0,1, 2, 3,.... 
h 


( 8 ) 


Now, solving the ordinary differential equation of X 
yields (noting the double root at m = 0) 


M co . f - m * x 

/ \ 


+ + h 

C o X + E c m e COS 

m7ry 

k 


m=l 

k n J 




M co , f 22 

? \ 


+ 

C o + E c rn e h cos 

mxy 

k 



m=l 

k 11 J 



( 9 ) 


where M coef , the number of coefficients used in the 
expansion, must be set a priori. From now on, the eigen- 
function terms are called modes, as commonly used in 
acoustic theory. The two separate roots have been dis- 
tinguished by a + or — superscript. For large values of 
negative x, the upstream velocity boundary condition 
requires 


aT 



( 10 ) 


Consider the entrance portion of the duct as shown 
in Fig. 4 for x < x in that is upstream to the duct 
obstruction. Employing separation of variables for the 
solution of Eq. (3), 


<f> & = X(x) Y(y) (4) 


thus, all values of c* must be set to zero and Eq. (9) 
becomes 


M, 


coef 


V x + 


+ E 

m=l 


ce h 
m 


cos 


mxy 


(ii) 


yields 

i£x = __i#y = j! (5) 

X 

The subscript a denotes the analytical solution in the 
entrance portion of the duct, as labeled in Fig. 4. 
Solving the ordinary differential equation for Y and 
applying the solid wall boundary conditions, 


Physically, the constant c Q in Eq. (11) represents a 
negative potential that will later be shown to be propor- 
tional to the size of the obstacle in the duct. The last 
term in Eq. (ll) represents damped higher order modes 
that blend the distorted potential around the obstacle 
into the uniform potential upstream. 

For convenience, the constant c will be pulled into 
the exponential, so that the summation begins with 
m = 0 such that 
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( 12 ) 


Interface Conditions 


M eoef _ ^ 

K = K x + E c rn e h COS 

m=0 



Also, because the damped exponentials can be very 
small at the beginning of the analytical regions, the 
separation constants are redefined. By introducing the 
identity in the brackets [], 



- 



M cotf 

mwx- n mxx in 

mxx 

f 

+ E c m 

^ h _ h 

e e 

e h cos 

mny 
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II 

o 



K n ) 


(13) 


Defining the constant A m as 


A 


m 


c 

m 


e 


h 


(14) 


At the interface 3 between the finite element domain 
and the analytical domains, both the potential and 
velocity are continuous. The local continuity of potential 

*ls- - *l s - = 0 I 18 ) 


can be expressed numerically by a collocation procedure 
(Lee and Cendes, 1987) or an integral weighting proce- 
dure (Baumeister and Krieder, 1992). The latter are 
used here with weight functions cos(m 7ry), so the con- 
tinuity of potential at the entrance interface S a is 
expressed by 


y= l 

f s w m . [*.-*]*- 0. 

“s, u 

y=° (19) 

W m< = cc^m’jry) 


such that 


(m* = 0,1,2, ,M coef equations) 


M mT(x-x in ) , 

CO«f 


K = + E a ; 


m=0 


m^ry 


(15) 


Finally, the characteristic length is set equal to h^ 
so that h = 1 and Eq. (15) simplifies to 


M 


coef 


K = + E A 


— mjr(x— x in ) 


cos(m?ry) 


(16) 


m=0 


Similarly, for the exit duct where x > x Qut , the 
analytical expression for the potential becomes 


h = + E B m e m,r(X x ° ut)cos ( m,r y) ^ 

m=0 


In a typical analytical solution, the separation 
constants A” and B^ are evaluated by some ortho- 
gonality condition in coordinate systems that matches 
the physical boundary. The rectangular coordinate 
system chosen herein obviously does not match the cir- 
cular boundary of the obstacle. However, the numerical 
grid system will transfer the necessary information such 
that A” and B* can be conveniently determined. 


There must be one equation for every unknown separa- 
tion constant. Thus, M coef + 1 weight equations are 
required. The superscript * designates the particular 
weight index to differentiate it from the eigenvalue 
index m used in Eq. (8). In contrast to FE weighted 
residual theory using local weight functions, the weight 
function here is global in nature acting over the whole 
y domain. 

Rather than take advantage of the orthogonality 
conditions in Eq. (19), a completely numerical approach 
has been adopted to determine equations for the A* in 
anticipation of future code complications. It suffices to 
apply a simple quadrature to obtain acceptable results 
when approximating Eq. (19). S a is divided into sub- 
intervals centered at points (x Is ,y Is ), which correspond 
to finite element boundary nodes introduced later in the 
paper. The nodes are evenly spaced on the boundary 
with the index I g . Once the number of modes M coef has 
been assumed (based on convergence of the numerical 
solution) the grid is set up so that the number of nodes 
on S a is N s > 6M coef to adequately resolve all the har- 
monic terms in Eq. (16). 

Applying the trapezoidal rule to the chosen nodes 
gives 
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E S [^(Wl.) " ^i]cos(mV yi jAy It = 0 . 

I.=l 


( 20 ) 


(m* = 0,1,2, M coef equations) 

where Ay is the distance between nodes except at the 
two end points. At the ends, Ay is half the distance 
between nodes. Thus, 


1 ~ + ) 

. 2\ * * p‘*/ 

A yi. = T7 : (21) 


M pts - 1 


S = Kronecker A 


By expressing <j> & in terms of the modal coefficients 
in Eq. (16), Eq. (20) can be written explicitly as 


M, 




cocf m ptl 

E A m Es co *(m»y I# )«>«(m**y Ii )Ay Ii 

m =° T=x 


M 


pt* 


-Es cos(m ,, 5ryj )<j > j Ayj 
1=1 


(22) 


M 


pt* 

"\ 

I =1 


“ 5^ cos(m*»ry 1> )Ay 1< 


(m*= 0,l,2,...,M coef equations) 


matrix system that yields values of all the unknowns 
values at the nodes as well as the separation constants 
AI and Bl. 


Finite Element Solution 


The finite element domain, with total area A, is 
divided into M discrete triangular elements, A^, 
e = defined by N corner nodal points (xr,yj), 

I = 1,2,...,N. The corner nodes for element area A' e ^ are 
denoted (x[ e ^,y[^), (x^,y^), and (x^*\y^). 

The potential is approximated by a linear combina- 
tion of weight functions Wj(x,y): 

?H x .y) = E w i( x >y)*i = [ w ( x >y)l W > ^ 

i=i 


with [ ] representing a row vector and { } representing 
a column vector. The weights have the property that 

Wi(xj, yj ) = *u (Kronecker tf) , ( 25 ) 


so that the unknown nodal pressure values are given by 

h = <H x i.yi)- 

To determine {^}, apply the method of weighted 
residuals. In this method, the residual error of Eq. (3), 


R r = J J A Wj (v • V* )dx dy 

(I = 1,2,3,,..,N one residual 


(26) 


Equation (22) comprises M coef + 1 separate difference 
equations, each of which is written in terms of all the 
unknown coefficients A m and the potential <^j s at the 
nodes on S a . 

The continuity of velocity requires that at the 
interface 


a* a 


dn 


Is* 



(23) 


equation for each node I) 

is set to 0 for each node I. Applying Green’s vector 
identity (integration by parts) and the divergence 
theorem to the integrand in Eq. (26) yields the weak 
formulation of Eq. (3) : 

Rj = J J A (vWj • V* )dx dy - JgfWjW ■ n)ds = 0 . 

= 1,2,3,...,N ) 



where n is the outward normal. This relationship is used 
to satisfy the natural boundary condition later required 
in the finite element portion of the problem. 

A similar set of equations can be written for the 
outlet region coefficients B^. These equations are later 
combined with the finite element equations to form a 


Equation (27) is a global, or node-oriented, formulation, 
in that it provides a difference equation for each node 
that can be used to determine {<j>}. 

From a practical standpoint, though, it is more 
convenient to consider a local, or element-oriented, 
formulation. To develop the local formulation, write 
each residual Rj as the sum of the element residuals: 


M M 

0 = R|=ERj*'= E 

e=l e=l 


JJawH” ■ ▼»«)** 0=R '-,I l^<-> K , '^" , ) dA 



V^ (e) • ii 
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(I as 1,2, 3, ...,N) f ( 28 ) 


where is the boundary of element A' e ). 


In the boundary integral terms in Eq. (28), it is rea- 
sonable to approximate the (continuous) normal deriva- 
tive with its mean value over S^nS. The key step is to 
apply the continuity of velocity (Eq. (23)) , which intro- 
duces the eigenfunction coefficients, thus linking the 
analytic solution and the finite element solution on the 
interface. The term is transformed as follows: 



(I = 1,2, 3,..., N) 


(31) 


The gradient of the potential in Eq. (31) can be 
evaluated directly from Eqs. (16) and (17) in terms of 
the unknown amplitude coefficients A and B + . 




ins 


W w^ (e) 
1 "aT 


Ids 




V *1/ 


's<*> 


f ( w ! c) ds 

Js (e) ns_ 1 


+ 


* b 


dn ) 


>S<«) 


f Mn W ( T e) ds, 

Jswns. i 


(29) 


Only the entrance and exit interfaces contribute, since 
the normal derivative of the potential is zero along the 
upper and lower channel walls and the obstacle. 

If fi is the angle between the positive x axis and the 
outward normal, for the geometry shown in Fig. 4, the 
interface between the elements and the analytical regions 
is vertical and fi = 180° at the entrance; thus, the 
normal derivative can be simplified to 


A 


dS 33, 

__cos(/3) + __sin(/3) 




,an J 

S<«) 

lax ay J 

s(«) 

. ax , 


(30) 

Similar for the exit except fi = 0 so that the cos(fi) has 
a value of +1. 

Substituting Eqs. (29) and (30) into Eq. (28) yields 


To evaluate the integrals in Eq. (31), it is necessary 
to represent ^(x,y) locally. Let Nj e ), j = 1,2,3, be the 
local linear shape functions for linear triangular elements 
associated with each corner node (Segerlind, 1976, 
p. 29), so that 

* (e) (x,y) = Nj'W)^ + N< e) (x,y) 

+ N< e) (x,y ) 4 e) 

= EN; t, (x,y)^ = |NW(x,y)]{^}. 
j = l 

(32) 

Next, implement the Galerkin method. When the global 
index I equals the local node index i associated with 
node (x[ e ),yj e )), let w[ c ) equal the local shape function 
N^. The global shape function Wj is assumed to be 
identically zero for any element where node I does not 
appear (simple pyramid weight approximation); thus, 
the line integral in Eq. (28) vanishes unless node I is on 
the boundary S. 

The solution to Eq. (31) is obtained by performing 
the usual element by element formation of the global 
matrix as presented in introductory FE texts. The Final 
form of the global matrix equation is obtained by com- 
bining the solution of Eq. (31) along with Eq. (22) and 
the exit equivalent to Eq. (22). 

[K] {$} = {F}, (”) 
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where 


(*} T = 


A o ,A 1 ,A 2 1 ‘"’ A M eo ,f’ 


(34) 


F contains the free stream velocity terms present in 
Eq. (22) and from the free stream velocity terms that 
enter Eq. (31) through the derivatives of <^ a and <£ a . 

The matrix K has the following general form 


A UU ! ^ 0 

A ( 21 ) : $( 22 ) i b< 23 > 

o : $( 32 ) | b< 33 ) 


This condition is necessary because complete information 
about the magnitude of the analytical duct modes was 
not passed to the nodal difference equations. Information 
about the entrance condition is passed to the finite ele- 
ment equations through the x derivative of the analyti- 
cal solution in Eq. (31). However, the magnitudes of the 
lowest order reflected mode A 0 and lowest order trans- 
mitted mode B^ are not passed into the modal differ- 
ence equations because the normal derivatives of the 
lowest order modes are identical to zero. 

In retrospect, Eq. (36) was not required in acoustics 
and electromagnetic applications of the modal element 
method discussed in the introduction. In these applica- 
tions, the Helmholtz equation governs rather than the 
Laplace equation. Consequently, the reflected and trans- 
mitted lowest order modes are harmonic wave like func- 
tions of x and derivatives exist for all modes. Thus, 
sufficient information is transmitted to the finite ele- 
ment equations to make them nonsingular. 


(35) Geometrical Model and Exact Solutions 


The top and bottom rows in the matrix represent the 
contribution from Eq. (22) for the entrance and its 
equivalent exit representation. The middle row repre- 
sents the contribution from Eq. (31). 

The submatrix A^ 11 ^ is a full M coe f X M coef matrix 
composed of the coefficients of the A m terms in Eq. (22) 
which resulted in applying Eq. (18) at the entrance 
interface. The submatrix $( 12 ^ is a sparse M coef x N 
matrix composed of the coefficients of <f > j in Eq. (22). 
The submatrices $^ 32 ^ and B^ 33 ^ are similar in form 
resulting from applying Eq. (18) at the exit interface. 

The submatrix A^ 21 ^ is an N x M coef matrix com- 
posed of the coefficients of A m from the surface integral 
in Eq. (31). For each boundary node, there is a full row 
of terms in A^ 21 ^, For interior nodes not on the bound- 
ary, there is a full row of zeros in A^ 21 ^ since the surface 
integral in Eq. (31) does not contribute to the difference 
equation at the interior nodes. A similar interpretation 
applies to B^ 23 ^. $( 22 ^ is a sparse, highly banded N x N 
matrix composed of the coefficients of <f > j resulting from 
the solution of Eq. (31). Equation (33) is solved by a 
standard frontal solver program. 

Dirichlet Condition 


One additional constraint is required to keep the 
matrix (Eq. (35)) nonsingular; namely the potential 
must be given a value (grounded) at some nodal location 
in the finite element grid. For the purposes of this paper, 

— 0.0 at x = 0.0 y = 1.0 (36) 


The cylindrical like obstacle shown in Fig. 2 is 
described by an infinite row of doublets in a uniform 
stream. Here, the complex potential can be written as 
(Kirchhoff, 1985) 


W(z) = ^(x,y) -b itf(x,y) = Ujs + 


*b 2 U + 


!coth 


'K 

12 ) 


(37) 


where 


z = x -b iy 


(38) 


and where the notation has been modified to reflect the 
nondimensionalization used herein. Using the identity 
given in Abramowitz and Stegun (1964, p. 84), the 
potential and stream function can be written as 


= U + x + 

~ CO 


* h2 K sinh(xx) 

2 cosh(^x) — cos(7ry) 


(39) 


* = U > - 


>rb 2 lf 

2 


sin(jry) 

cosh(jrx) — cos(^y) 


(40) 


In the examples to follow, 

U + = 1.0 ( 41 ) 
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so that the value of the streamline along the upper wall 
y — 1 is 0 = 1 and the value of the streamline is t/> = 0 
along the lower wall and obstacle. Therefore, the height 
of the obstacle y as a function of x is given by 


xb 2 sin(jry) 

2 cosh(7rx) - cos(fl-y) 


(42) 


sinh(7rx) 

cosh(^rx) — cos(xy) 

OB 

= -1 - 2 e m * x cos(m»ry) (44) 
m=l 


X < 0 


Equation (42) is solved numerically by the bisection 
method for the height y of the obstacle as a function x 
and as a function of the assumed strength of the dipole 
parameter b. The following two values of b are used in 
the examples of the next section: 

b x o y 0 

0.5042 0.5084 0.5000 

1.9020 1.1552 0.9000 

Recall that the extremes of the obstacle x 0 and y 0 are 
defined in Fig. 2. For b equals 0.5642 the obstacle 
closely approximates a circular cylinder. 

Equation (39) represents the exact analytical solution 
to which the approximate numerical analysis will be 
compared graphically in the examples to follow. In addi- 
tion, the numerically determined modal coefficients in 
Eq. (16) can also be compared to the exact results in 
Eq. (39) by equating both expressions: 

A 0+ E A m e cos(mjry) 

m=1 (43) 

__ jrb 2 sinh(jrx) 

2 cosh(7rx) — cos(?ry) 

The lowest order reflected mode A 0 can be easily deter- 
mined by letting x go to negative infinity. Noting that 
the higher order modes on the left-hand side go to zero 
and the ratio of sinh to cosh is — 1 at — it follows 



Thus, the lowest order potential acts with the opposite 
polarity of the free stream potential and is proportional 
to the square of the dipole strength. 

The amplitude of the higher order modes can also be 
determined using the series expansion (Gradshteyn and 
Ryzhik, 1965, p. 42, Eq. (1.461), modified for x < 0) 


Substituting Eq. (44) into Eq. (43) , the exact expres- 
sion for the modal amplitude A m can be expressed as 

A" = - J rb 2 e mTXi “ («) 

m 

Equations (43) and (45) can be used to compare the 
exact results with the analytical approximations. A simi- 
lar analysis applies to the B+ modes. For this symmetri- 
cal obstacle the modes are identical in magnitude to 
the A~ modes but of opposite sign. 

For the examples to be considered later, the follow- 
ing tabulated values give the exact amplitude of the 
back reflected potential A 0 and the higher order cutoff 


potential modes. 


Example 1 

Example 2 

b = 0.5042 

b = 1.9020 

x 0 = -0.5084 

x = -1.1552 

o 

x in = -0.5084 

x in = -1.3984 

Aq = — 0.5000E+00 

A^ = — 5.0825E+OO 

A7 = — 0.2024E+00 

AY = -0.1404E+00 

A^ = — 0.4099E— 01 

AY = — 0.1736E— 02 

A~ = — 0.8299E— 02 

AY = —0.2146E— 04 

A7 = — 0.1680E— 02 

AY = — O.2052E— 06 

A~ = — 0.3401E— 03 

AY = — 0.3279E— 08 

A^ = — O.0887E— 04 

AY = — 0.4053E— 10 


As seen above, the amplitudes of the higher order 
modes fall off rapidly with increasing m. Thus, the grid 
density can be reasonably sparse in the transverse y 
direction to resolve the important modes. Because the 
grid was extended farther from the body in the second 
example (x in < x Q ), the modal element coefficients for 
the higher order modes are smaller for this second exam- 
ple. Of course, as the magnitude of x in increases to very 
large values, all the higher order modes will become neg- 
ligible, so the grid density could be very sparse near the 
end of the finite element grid in the neighborhood x in . 
However, the analysis also indicates that the finite ele- 
ment grid density must still be increased at x Q to what- 
ever level of accuracy is desired to resolve the higher 
order modes. 
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Results and Comparisons 

To validate the method, two numerical experiments 
are presented for potential flow over a cylinder like 
obstacle described by Eq. (42). The first example con- 
siders an obstacle that blocks 50 percent of the channel 
while the second example is concerned with an obstacle 
blocking 90 percent of the channel. For the examples to 
follow, recall that the characteristic length is set 
equal to the height of the channel so that the channel 
has a dimensionless height of unity. The dimensionless 
mean flow velocity is also assumed to be unity. The 
parameter M coef specifying the number of modes in 
Eq. (8) was taken to be 3 for a good comparison to the 
exact solution. For problems without an exact solution, 
the number of modes must be increased till the results 
converge. For example, in the acoustic application 
shown in Fig. 1(d), the number of modes required was 
121 . 

Example 1. — Half Channel Obstruction 


potential were calculated in each analytical region as 
shown in Fig. 5. Clearly, the modal element method 
gives good agreement with the numerical results. 

A convergence check was made in this example by 
increasing the number of vertical nodes to 25 and x 
coordinate nodes to 21 for a total of 525 nodes and 960 
elements. In this case, the numerically calculated modes 
and the exact modes are 

Numerical Solution Exact Solution 


A“ = — 0.497E— 00 A“ = -0.50002E-00 

AY = — 0.202E— 00 AY = -0.20247E-00 

AY = — 0.446E— 01 AY = -0.40994E-01 

No improvement of the graphical results was seen by 
eye. Figure 6 shows the resulting contour plot including 
the analytical and finite element regions. The dash line 
in Fig. 6 shows the streamlines. There is good agreement 
between the exact and the modal element results. 


Consider the potential flow over the cylinder with a 
b value of 0.5642, x 0 = —0.5082, and y 0 = 0.5000. In 
this example, the finite element region was placed 
directly over the obstacle as shown in Fig. 4 so that the 
end of the analytical domain x in coincides with the 
beginning of the obstacle at x Q . As shown in Fig. 4, the 
finite element domain extends from -0.5082 to 0.5082. 
Thirteen nodes were used on the interface and eleven 
nodes along the surface of the obstacle for a total of 
143 nodes and 240 elements. 


In Fig. 5, the velocity potential is plotted as a 
function of x along the upper wall, y = 1. The dashed 
line represents the potential on the duct wall without an 
obstacle. The modal element solutions (hollow boxes) 
are compared to the exact solution (solid line) given by 
Eq. (39). In Fig. 5, for -0.5082 < x < 0.5082, the 
values of potential at the finite element nodes are used 
to generate the solution. Here, eleven closely packed 
nodal values are shown in Fig. 5. 

Also in Fig. 5, the numerical solution is generated 
from Eq. (16) using the numerically determined modal 
coefficients A^ for x < —0.5082 and for x > 0.5082. 
The numerically calculated modal coefficients and the 
exact coefficients are: 


Numerical Solution 

AY = — 0.486E— 00 
AY = —0.197E— 00 
AY = — 0.454E— 01 


Exact Solution 

AY = — 0.50002E— 00 
AY = — 0.20247E— 00 
AY = — 0.40994E— 01 


The disagreement in the AY modal amplitude was 
believed to be a result of the highly skewed triangles 
near the leading edge of the cylinder as shown in Fig. 3. 
This simple grid system was chosen just to confine the 
grid directly over the obstacle. In the next example with 
90 percent channel blockage and a very steep slope near 
the leading edge, a more conventional grid system is 
employed. As will now be shown, the new grid system 
will lead to good resolution of the highest order mode. 

Example 2. — Large Channel Obstruction 

Consider the potential flow over the cylinder with a 
b value of 1.9020, x 0 = 1.1552, and y G = 0.9000. In this 
case, 90 percent of the channel has an obstruction as 
illustrated in Fig. 7. The grid was extended slightly in 
front of the obstacle to better resolve the steep slope of 
the obstacle near x 0 . In this case, x in = —1.3984 and 
x out = 1.3984. Twenty one nodes were used on the inter- 
face and 132 nodes along the upper channel wall for a 
total of 1272 nodes and 2212 elements, as shown in 
Fig. 7. 

In Fig. 8, the velocity potential is again plotted as a 
function of x along the upper wall, y = 1. The dashed 
line again represents the potential on the duct wall with- 
out an obstacle. The modal element solutions (hollow 
boxes) are compared to the exact solution (solid line) 
given by Eq. (39). In Fig. 8, for —1.3984 < x < 1.3984, 
the values of potential at the finite element nodes are 
used to generate the solution. In this case, the closely 
packed nodal values are shown in Fig. 8. 


The numerically calculated and exact modal coefficients Also in Fig. 8, the numerical solution is generated 

are in reasonable agreement. Six separate values of the from Eq. (16) using the numerically determined modal 
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coefficients A~ for x < —1.3984 and B* for x > 1.3984. 
The numerically calculated modal coefficients and the 
exact coefficients are as follows: 


Numerical Solution Exact Solution 


A” = — 0.567E+01 A o = -0.56825E+01 

A 7 = —0. HOE— 00 A 7 = — 0.14048E— 00 

A7 = — 0.173E— 02 A 7 = -0.17364E-02 

In this example, the numerically calculated and exact 
modal coefficients are in good agreement for all modes 
including the highest order mode. Six separate values of 
the potential were calculated in each analytical region as 
shown in Fig. 8. Clearly, the modal element method 
again gives good agreement with the numerical results. 


Finally, Fig. 9 shows a contour plot of the potential 
inside the finite element region while Fig. 10 shows a 
contour plot including both the analytical and finite ele- 
ment regions. The dash line in Fig. 10 shows the stream- 
lines. Again, there is good agreement between the exact 
and the modal element results. 


exact magnitude of the back potential and the decay 
rates of the harmonics which blend the flow streamline 
and the potential lines from about the obstacle into the 
uniform flow lines of the far field. 

Eigenfunction solutions are applicable for a wide 
range of practical CFD problems in regions where viscos- 
ity no longer dominates. Nevertheless, eigenfunction 
solutions do not exist for most CFD problems. For these 
more complicated problems, a challenging and intriguing 
aspect of the modal element approach could be to use a 
finite series of known trial functions with unknown coef- 
ficients Aj. These trial functions would approximate the 
physics, but not necessarily satisfy the governing differ- 
ential equations. Meirovitch (1967) suggests using 
admissible trial functions (satisfying geometric boundary 
conditions) or comparison trial functions (satisfying 
geometric and natural boundary conditions) for this 
task. In these cases, additional constraints must be 
applied to the coefficients Aj so that differential equa- 
tions are satisfied in the analytical region. Meirovitch 
suggests that the collocation method is a practical 
approach to this problem. 


The first mode A 7 has increased from —0.5 to 
—5.6825 or a factor of 10 increase in the magnitude of 
the potential. This is a direct result of the larger 
obstruction in the second example. In both cases the 
gradient of the potential along the upper wall reaches 
the free stream value once outside the obstacle because 
of the quick decay of the higher order modes. 

Concluding Remarks 

The modal element method for potential flow over a 
two-dimensional cylindrical like obstacle is presented. 
The total flow domain is broken into three subdomains 
that are patched together. The potential field is repre- 
sented by a finite element solution in the irregular 
subdomain next to the obstacle and by an exact eigen- 
function expansion in the unbounded entrance and exit 
ducts. The analytical and numerical solutions are 
coupled by the continuity of potential and velocity 
across the interface between the subdomains and are 
calculated simultaneously from a single matrix equation. 
The method is applicable to problems involving a com- 
plete range of channel blockage. 

The combined numerical and analytical results show 
excellent agreement with the corresponding exact solu- 
tions. For numerical insight, the analytical results indi- 
cate the accuracy of the chosen exit boundary condition 
and the grid density required for a given harmonic 
accuracy near the obstacle (generally about 12 nodes per 
wavelength are required to resolve a modal harmonic). 
For flow field insight, the analytical results indicate the 


The long term goal or vision of this research is to 
(1) adapt the modal element method for greater analyti- 
cal and numerical insight to a wider class of CFD prob- 
lems and (2) decrease computational costs by reducing 
the numerical grid. 
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(a) Computational grid; Khan, Brown, Ahuja 
transient difference. 



radius of cylinder = 5 cm 


(b) Computational grid; model ring grid analysis. 


Exact analysis ABS (p) 

• Khan, Brown, Ahuja transient difference 

□ Modal ring grid analysis 



(c) Comparison of Numerical Approaches k# = 
0.182 cm' 1 radius of cylinders = 5 cm, k*r" * 1 .08. 


Figure 1 . — Application of modal element methods 


Exact analysis 

□ Numerical solution 



in acoustic scattering from a hard circular cylinder. 
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Figure 3. — Conventional finite element discretization for flow around 
cylindrical object in channel. 
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and numerical 
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Figure 4.— Modal element finite element discretization for flow around cylindrical object 
in channel. 
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x axial coordinate 


Figure 5.— Effect of a half channel obstruction on the potential 
along the upper wall (143 nodes and 240 elements). 


Exact analysis 

□» etc Numerical solution 



Figure 6. — Contour plots of the potential in both the finite element region and the 
analytical regions for the half channel obstruction (525 nodes and 960 elements). 



X 


out 


Figure 7. — Modal element finite element discretization for flow around large cylindrical object 
in channel. 
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x axial coordinate 

Figure 8. — Effect of a large channel obstruction on the potential 
along the upper wall (1 272 nodes and 221 2 elements). 



- 2-10 1 2 


Shape of boundary 

Figure 9. — Contour plots of the potential in the finite element 
region for the channel obstruction. 


Exact analysis 

□ , etc Numerical solution 



Figure 10. — Contour plots of the potential In both the finite element region 
and the analytical regions for the half channel obstruction (1 272 nodes and 
221 2 elements). 
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